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Abstract 

This paper describes the algorithms, features and implementation of PyDEC, a Python library for 
computations related to the discretization of exterior calculus. PyDEC facilitates inquiry into both 
physical problems on manifolds as well as purely topological problems on abstract complexes. We 
describe efficient algorithms for constructing the operators and objects that arise in discrete exte- 
rior calculus, lowest order finite element exterior calculus and in related topological problems. Our 
algorithms are formulated in terms of high-level matrix operations which extend to arbitrary dimen- 
sion. As a result, our implementations map well to the facilities of numerical libraries such as NumPy 
and SciPy. The availability of such libraries makes Python suitable for prototyping numerical meth- 
ods. We demonstrate how PyDEC is used to solve physical and topological problems through several 
concise examples. 

Categories and Subject Descriptors: G.4 [Mathematical Software] ; G.l .8 [Numerical Analysis] : 
Partial Differential Equations - Finite element methods, Finite volume methods, Discrete exterior cal- 
culus, Finite element exterior calculus; 1.3.5 [Computer Graphics] : Computational Geometry and Ob- 
ject IVlodeling- Geometric algorithms, languages, and systems, Computational topology 

1 Introduction 

Geometry and topology play an increasing role in the modern language used to describe physical prob- 
lems [ 1 , 25] . A large part of this language is exterior calculus which generalizes vector calculus to smooth 
manifolds of arbitrary dimensions. The main objects are differential forms (which are anti-symmetric 
tensor fields), general tensor fields, and vector fields defined on a manifold. In addition to physical ap- 
plications, differential forms are also used in cohomology theory in topology [14] . 

Once the domain of interest is discretized, it may not be smooth and so the objects and operators of 
exterior calculus have to be reinterpreted in this context. For example, a surface in ^ may be discretized 
as a two dimensional simplicial complex embedded in K^, i.e., as a triangle mesh. Even when the domain 
is a simple domain in space, such as an open subset of the plane or space the discretization is usually 
in the form of some mesh. The various objects of the problem then become defined in a piecewise 
varying fashion over such a mesh and so a discrete calculus is required there as well. After discretizing the 
domains, objects, and operators, one can compute numerical solutions of partial differential equations 
(PDEs), and compute some topological invariants using the same discretizations. Both these classes of 
applications are considered here. 

There have been several recent developments in the discretization of exterior calculus and in the clar- 
ification of the role of algebraic topology in computations. These go by various names, such as covolume 
methods [46], support operator methods [50], mimetic discretization [10, 36], discrete exterior calculus 



(DEC) [19, 27, 31], compatible discretization, finite element exterior calculus [2, 4, 35], edge and face 
elements or Whitney forms [13, 30], and so on. PyDEC provides an implementation of discrete exterior 
calculus and lowest order finite element exterior calculus using Whitney forms. 

Within pure mathematics itself, ideas for discretizing exterior calculus have a long history. For exam- 
ple the de Rham map that is commonly used for discretizing differential forms goes back at least to [16] . 
The reverse operation of interpolating discrete differential forms via the Whitney map appears in [52] . A 
combinatorial (discrete) Hodge decomposition theorem was proved in [21] and the idea of a combina- 
torial Hodge decomposition dates to [23] . JVEore recent work on discretization of Hodge star and wedge 
product is in [53, 54] . Discretizations on other types of complexes have been developed as well [28, 49] . 

1 . 1 Main contributions 

In this paper we describe the algorithms and design of PyDEC, a Python software library implementing 
various complexes and operators for discretization of exterior calculus, and the algorithms and data 
structures for those. In PyDEC all the discrete operators are implemented as sparse matrices and we 
often reduce algorithms to a sequence of standard high-level operations, such as sparse matrix-matrix 
multiplication [6], as opposed to more specialized techniques and ad hoc data structures. Since these 
high-level operations are ultimately carried out by efficient, natively-compiled routines (e.g. C or Fortran 
implementations) the need for further algorithmic optimization is generally unnecessary. 

As is commonly done, in PyDEC we implement discrete differential forms as real valued cochains 
which will be defined in Section 2. PyDEC has been used in a thesis [8], in classes taught at University of 
Illinois, in experimental parts of some computational topology papers [20, 22, 32], inDarcy fiow [34], and 
in least squares ranking on graphs [33]. The PyDEC source code and examples are publicly available [7]. 
We summarize here our contributions grouped into four areas. 

Basic objects and metbods: (1) Data structures for : simplicial complexes of dimension n embedded 
in R^, N'>n; abstract simplicial complexes; Vietoris-Rips complexes for points in any dimension; and 
regular cube complexes of dimension n embedded in IR"; (2) Cochain objects for the above complexes; 
(3) Discrete exterior derivative as a coboundary operator, implemented as a method for cochains on 
various complexes. 

Finite element exterior calculus: (1) Fast algorithm to construct sparse mass matrices for Whitney forms 
by eliminating repeated computations; and (2) Assembly of stiffness matrices for Whitney forms from 
mass matrices by using products of boundary and stiffness matrices. Note that only the lowest order 
(S^j") elements of finite element exterior calculus are implemented in PyDEC. 

Discrete exterior calculus: (1) Diagonal sparse matrix discrete Hodge star for well-centered (circum- 
centers inside simplices) and Delaunay simplicial complexes (with an additional boundary condition); 
(2) Circumcenter calculation for fc-simplex in an ^-dimensional simplicial complex embedded in 
using a linear system in barycentric coordinates; and (3) Volume calculations for primal simplices and 

circumcentric dual cells. 

Examples: (1) Resonant cavity curl-curl problem; (2) Flow in porous medium modeled as Darcy flow, i.e., 
Poisson's equation in first order (mixed) form; (3) Cohomology basis calculation for a simplicial mesh, 
using harmonic cochain computation using Hodge decomposition; (4) Finding sensor network coverage 
holes by modeling an abstract, idealized sensor network as a Rips complex; and (5) Least squares ranking 
on graphs using Hodge decomposition of partial pairwise comparison data. 
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2 Overview of PyDEC 



One common type of discrete domain used in scientific computing is triangle or tetrahedral mesh. These 
and their higher dimensional analogues are implemented as n-dimensional simplicial complexes em- 
bedded in IR^, N> n. Simplicial complexes are useful even without an embedding and even when they 
don't represent a manifold, for example in topology and ranking problems. Such abstract simplicial 
complexes without any embedding for vertices are also implemented in PyDEC. The other complexes 
implemented are regular cube complexes and Rips complexes. Regular cubical meshes are useful since 
it is easy to construct domains even in high dimensions whereas simplicial meshing is hard enough in 
3 dimensions and rarely done in 4 or larger dimensions. Rips complexes are useful in applications such 
as topological calculations of sensor network coverage analysis [17]. The representations used for these 
four types of complexes are described in Section 3-6. A complex that is a manifold (i.e., locally Euclidean) 
will be referred to as mesh. 

The definitions here are given for simplicial complexes and generalize to the other types of complexes 
implemented in PyDEC. In PyDEC we only consider integer valued chains and real-valued cochains. 
Also, we are only interested in finite complexes, that is, ones with a finite number of cells. Let K be a 
finite simplicial complex and denote its underlying space by \K\. Give \K\ the subspace topology as a 
subspace of IR^ (a set U m\K\ is open iff Un\K\is open in IR^). For a finite complex this is the same as 
the standard way of defining topology for |^| [45, pages 8-9] and \K\isa closed subspace of M^. 

An oriented simplex with vertices vo,...,Vp will be written as [voy-> Vp] and given names like 
with the superscript denoting the dimension and subscript denoting its place in some ordering of p- 
simplices. Sometimes the dimensional superscript and/or the indexing subscript will be dropped. The 
orientation of a simplex is one of two equivalence classes of vertex orderings. Two orderings are equiv- 
alent if one is an even permutation of the other. For example [i^o. v\, V2] and [1^1, V2, vq] denote the same 
oriented triangle while [1^0, ^2. ^"1] is the oppositely oriented one. 

A p-chain of K is a function c from oriented p-simplices of K to the set of integers Z, such that 
c{-a) = -cia] where -a is the simplex a oriented in the opposite way. Two chains are added by adding 
their values. Thus p-chains are formal linear combinations (with integer coefficients) of oriented p- 
dimensional simplices. The space of p-chains is denoted Cp[K) and it is a fi:ee abelian group. See [45, 
page 21]. Free abelian groups have a basis and one does not need to impose a vector space structure. For 
example, a basis for CpiK) is the set of integer valued functions that are 1 on a p-simplex and on the 
rest, with one such basis element corresponding to each p-simplex. These are called elementary chains 
and the one corresponding to a p-simplex will also be referred to as af. The existence of this basis 
and the addition and negation of chains is the only aspect that is important for this paper. The intuitive 
way to think of chains is that they play a role similar to that played by the domains of integration in the 
smooth theory. The negative sign allows one to talk about orientation reversal and the integer coefficient 
allows one to say how many times integration is to be done on that domain. 

Sometimes we will need to refer to a dual mesh which will in general be a cell complex obtained from 
a subdivision of the given complex K. We'll refer to the dual complex as ★ K. For a discrete Hodge star 
diagonal matrix of DEC, the dual mesh is the one obtained from circumcentric subdivision of a well- 
centered or Delaunay simplicial complex and such a Hodge star is described in Section 10. 

Homomorphisms from the p-chain group Cp{K) to IRare called p -cochains of K and denoted C^{K; IR). 
This set is an abelian group and also a vector space over IR. Similarly the dual p-cochains are denoted 
C''(*^r;IR) or DP{J^K;U). The discretization map from space of smooth p-forms to p-cochains is called 
the rfei?/zam map R:QP(ir) ^ CP{K;m orR:nP{K) — CPi-kK-M)- See [16,21]. For a smooth p-form a, 
the de Rham map is defined as R : a (c f^a) for any chain c e CpiK). We will denote the evaluation 
of the cochain R(a) on a chain c as <R(a), c>. A basis for Cf{K;U) is the set of elementary cochains. The 
elementary cochain (cr^)* is the one that takes value 1 on elementary chain cr'' and on the other ele- 
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mentary chains. Thus the vector space dimension of CP{K;R] is the number of p-simphces in K. We'll 
denote this number by Np. Thus Nq will be the number of vertices, Ni the number of edges, N2 the 
number of triangles and so on. 

Like most of the numerical analysis literature mentioned in Section 1 we assume that the smooth 
forms are either defined in the embedding space of the simplicial complex, or on the complex itself, or 
can be obtained by pullback from the manifold that the complex approximates. In contrast, most math- 
ematics literature quoted including [21, 52] uses simplicial complex defined on the smooth manifold as 
a "curvilinear" triangulation. In the applied literature, the complex approximates the manifold. Many fi- 
nite element papers deal with open subsets of the plane or IR^ so they are working with triangulations of a 
manifold with piecewise smooth boundaries. Surface finite element methods have been studied outside 
of exterior calculus [18]. A variational crimes methodology is used for finite element exterior calculus 
on simplicial approximations of manifolds in [35] . In the computer graphics literature, piecewise-linear 
triangle mesh surfaces embedded in are common and convergence questions for operators on such 
surfaces have been studied [29] . In light of all of these, PyDEC's framework of using simplicial or other 
approximations of manifolds is appropriate. 

Operators such as the discrete exterior derivative (d) and Hodge star (*) can be implemented as 
sparse matrices. At each dimension, the exterior derivative can be easily determined by the incidence 
structure of the given simplicial mesh. For DEC the Hodge star is a diagonal matrix whose entries are de- 
termined by the ratios of primal and dual volumes. Care is needed for dual volume calculation when the 
mesh is not weU- centered. For finite element exterior calculus we implement Whitney forms. The cor- 
responding Hodge star is the mass matrix which is sparse but not diagonal. One of the stiffness matrices 
can be obtained from it by combining it with the exterior derivative. 

Once the matrices implementing the basic operators have been determined, they can be composed 
together to obtain other operators such as the codifferential (5) and Laplace-deRham (A). While this 
composition could be performed manually, i.e. the appropriate set of matrices combined to form the 
desired operation, it is prone to error. In PyDEC this composition is handled automatically. For exam- 
ple, the function d(.) which implements the exterior derivative, looks at the dimension of its argument to 
determine the appropriate matrix to apply. The same method can be applied to the codifferential func- 
tion which then makes their composition 5(d(.)) work automatically. This automation eliminates a 
common source of error and makes explicit which operators are being used throughout the program. 

PyDEC is intended to be fast, flexible, and robust. As an interpreted language, Python by itself is 
not well-suited for high-performance computing. However, combined with numerical libraries such as 
NumPy and SciPy one can achieve fast execution with only a small amount of overhead. The NumPy 
array structure, used extensively throughout SciPy and PyDEC, provides an economical way of storing 
N-dimensional arrays (comparable to C or Fortran) and exposes a C API for interfacing Python with other, 
potentially lower-level libraries [5 1] . In this way. Python can be used to efficientiy "glue" different highly- 
optimized libraries together with greater ease than a purely C, C-n-, or Fortran implementation would 
permit [47]. Indeed, PyDEC also makes extensive use of the sparse module in SciPy which relies on 
natively-compiled C++ routines for all performance-sensitive operations, such as sparse matrix-vector 
and matrix-matrix multiplication. PyDEC is therefore scalable to large data sets and capable of solving 
problems with millions of elements [8] . 

Even large-scale, high-performance libraries such as Trilinos provide Python bindings showing that 
Python is useful beyond the prototyping stage. We also make extensive use of Python's built-in unit 
testing framework to ensure PyDEC's robustness. For each non-trivial component of PyDEC, a number 
of examples with known results are used to check for consistency. 
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2. 1 Previous work 



Discrete differential forms now appear in several finite element packages such as FEMSTER [15] , DOLFIN 
[43] and deal.II [5] . These libraries support arbitrary order conforming finite element spaces in two and 
three dimensions. In contrast, for finite elements PyDEC supports simplicial and cubical meshes of ar- 
bitrary dimension, albeit with lowest order elements. In addition, PyDEC also supports the operators of 
discrete exterior calculus and complexes needed in topology. We note that Exterior [41], an experimen- 
tal library within the FEniCS [42] project, realizes the framework developed by Arnold et al. [3] which 
generalizes to arbitrary order and dimension. Exterior uses symbolic methods and supports integration 
of forms on the standard simplex. PyDEC supports mass and stiffness matrices on simplicial and cubi- 
cal complexes. The discovery of lower dimensional faces in a complex and the computation of all the 
boundary matrices is also implemented in PyDEC. 

The other domain where PyDEC is useful is in computational topology. There are several packages in 
this domain as well, and again PyDEC has a different set of features and aims from these. In [39] efficient 
techniques are developed for finding meaningful topological structures in cubical complexes, such as 
digital images. In addition to simplicial and cubical manifolds, PyDEC also provides support for abstract 
simplicial complexes such as the Rips complex of a point set. The Applied and Computational Topology 
group at Stanford University has been the source for several packages for computational topology. These 
include various versions of PLEX such as JPlex and javaPlex which are designed for persistent homol- 
ogy calculations. Another package from the group is Dionysus, a C++ library that implements persistent 
homology and cohomology [24, 55] and other interesting topological and geometric algorithms. In con- 
trast, we view the role of PyDEC in computational topology as providing a tool to specify and represent 
different types of complexes, compute their boundary matrices, and compute cohomology representa- 
tives with or without geometric information. 

3 Simplicial Complex Representation 

Before detailing the algorithms used to implement discretizations of exterior calculus, we discuss the 
representation of various complexes, starting in this section with simplicial complexes. Consider the 
triangle mesh shown in Figure 1 with vertices and faces enumerated as shown. This example mesh is 
represented by arrays 



where the subscript 2 denotes the dimension of the simplices. The i-th row of V contains the spatial 
coordinates of the i-th vertex. Likewise the i-th row of simplex array §2 contains the indicies of the 
vertices that form the i-th triangle. The indices of each simplex in §2 in this example are ordered in a 
manner that implies a counter-clockwise orientation for each. For an n-dimensional discrete manifold, 
or mesh, arrays V and S„ suffice to describe the computational domain. 

In addition to V and §„, an w-dimensional simplicial complex is comprised by its p-dimensional 
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Figure 1: Simplicial mesh with enumerated vertices and simplices. 
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Figure 2: Simplicial complex with oriented edges and triangles. 



faces, Sp,0 < p<n. In the case of Figure 1, these are 



So = 



which correspond to the vertices (0-simplices) and oriented edges (1 -simplices) of the complex. A graph- 
ical representation of this simplicial complex is shovm in Figure 2. Since the orientation of the lower (< n) 
dimensional faces is arbitrary, we use the convention that face indices will be in sorted order. Further- 
more, we require the rows of § to be sorted in lexicographical order. As pointed out in Section 7 and 9, 
these conventions facilitate efficient construction of differential operators and stiffness matrices. 



4 Regular Cube Complex Representation 

PyDEC provides a regular cube complex of dimension n embedded in U" for any n. As mentioned earlier, 
in dimension higher than 3, constructing simplicial manifold complexes is hard. In fact, even construc- 
tion of good tetrahedral meshes is still an active area in computational geometry. This is one reason 
for using regular cube complexes in high dimensions. Moreover, for some applications, like topological 
image analysis or analysis of voxel data, the regular cube complex is a very convenient framework [39] . 

A regular cube complex can be easily specified by an n-dimensional array of binary values (bitmap) 
and a regular n-dimensional cube is placed where the bit is on. For example the cube complex shown in 
Figure 3 can be created by specifying the bitmap array 

1 

1 1 ■ 
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Figure 3: Regular cube mesh with enumerated vertices and faces. 
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Figure 4: Regular cube complex with oriented edges and faces. 



A bitmap suffices to describe the top level cubes, but a cube array (like simplex array) is used during 
construction of differential operators and for computing faces. In tliis paper we describe the construc- 
tion of exterior derivative, Hodge star and Whitney forms on simplicial complexes. For cube complexes 
we describe only the construction of exterior derivative and lower dimensional faces. However, the other 
operators and objects are also implemented in PyDEC for such complexes. For example, Whitney-like 
elements for hexahedral grids are described in [11] and are implemented in PyDEC. 

Converting a bitmap representation of a mesh into a cube array representation is straightforward. 
For example, the cube array representation of the mesh in Figure 3 is 



C2 = 



1 

10 1 

11 1 



As with the simplex arrays, the rows of C2 correspond to individual two-dimensional cubes in the mesh. 
The two left-most columns of C2 encode the origins of each two-dimensional cube, namely (0,0), (1,0) 
and (1, 1). The remaining two columns encode the coordinate directions spanned by the cube. Since C2 
represents the top-level elements, all cubes span both the x (coordinate 0) and y (coordinate 1) dimen- 
sions. In general, the first n columns of Cjt encode the origin or corner of a cube while the remaining k 
columns identify the coordinate directions swept out by the cube. We note that the cube array represen- 
tation is similar to the cubical notation used by Sen [49] . 
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The edges of the mesh in Figure 3 are represented by the cube array 



Ci = 
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where again the first two columns encode the origin of each edge and the last column indicates whether 
the edge points in the x or y direction. For example, the row [0, 0, 0] corresponds to edge in Figure 4 
which begins at (0,0) and extends one unit in the x direction. Similarly the row [2, 1, 1] encodes an edge 
starting at (2, 1) extending one unit in the y direction. Since zero-dimensional cubes (points) have no 
spatial extent their cube array representation 



Co = 
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contains only their coordinate locations. 

The cube array provides a convenient representation for regular cube complexes. While a bitmap 

representation of the top-level cubes is generally more compact, the cube array representation general- 
izes naturally to lower-dimensional faces and is straightforward to manipulate. 



5 Rips Complex Representation 

The Rips complex, or Vietoris-Rips complex of a point set is defined by forming a simplex for every subset 
of vertices with diameter less than or equal to a given distance r. For example, if pair of vertices {Vi, vj) 
are no more than distance r apart, then the Rips complex contains an edge (1 -simplex) between the 
vertices. In general, a set of jo > 2 vertices forms a (p- 1) -simplex when all pairs of vertices in the set are 
separated by at most r. 

In recent work, certain sensor network coverage problems have been shown to reduce to finding 
topological properties of the network's Rips complex, at least for an abstract model of sensor networks [17] 
Such coordinate-free methods rely only on pairwise communication between nodes and do not require 
the use of positioning devices. These traits are especially important in the context of ad-hoc wireless 
networks with limited per-node resources. Figure 5 depicts a planar sensor network and its associated 
Rips complex. 

In this section we describe an efficient method for computing the Rips complex for a set of points. 
Although we consider only the case of points embedded in Euclidean space, our methodology applies to 
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Figure 5: Broadcast radii and Rips complex for a sensor network. 

more general metric spaces. Indeed, only the construction of the 1 -skeleton of the Rips complex requires 
metric information. The higher- dimensional simplices are constructed directly from the 1 -skeleton. 

We compute the 1-skeleton of the Rips complex with a kD-Tree data structure. Specifically, for each 
vertex vi we compute the set of neighboring vertices { Vj : || Vj - vi || < r}. The hierarchical structure of the 
kD-Tree allows such queries to be computed efficiently. 

The 1 -skeleton of the Rips complex is stored in an array §i , using the convention discussed in Section 
3. Additionally, the (oriented) edges of the 1-skeleton are used to define E, a directed graph stored in a 
sparse matrix format. Specifically, E{z, j) = 1 if [i, j] is an edge of the Rips complex, and zero otherwise. 
For the Rips complex depicted in Figure 6, 



§1 



1 

2 

3 

1 3 

2 3 



111 

1 

1 





Let Fp denote 



are the corresponding simplex array and directed graph respectively. 

The arrays of higher dimensional simplices §2.§3.-- - can be computed as follows 
the (sparse) matrix whose rows are identified with the p-simplices as specified by §p. Each row of Fp 
encodes the vertices which form the corresponding simplex. Specifically, fpii,j) takes the value 1 if the 
i-th simplex contains vertex j and zero otherwise. For the example shown in Figure 6, 



Fi 



110 
10 10 
10 1 
10 1 
11 



encodes the edges stored in §i. Once Fp is constructed we compute the sparse matrix-matrix product 
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Figure 6: Five directed edges form the 1 -skeleton of the Rips complex. 



Fp E. For our example the result is 



FiE = 



112 

112 

111 

1 

1 



Like Fp, the product Fp E is a matrix that relates the p-simplices to the vertices: the matrix entry [i, j) of 
Fp E counts the number of directed edges that exist from the vertices of simplex / to vertex j. When the 
value of {i, j) entry of Fp E is equal to p + 1, we form a p-simplex of the Rips complex by concatenating 
simplex / with vertex j. In the example, matrix entries (0, 3) and (1,3) of Fi E are equal to 2 which implies 
that the 2-skeleton of the Rips complex contains two simplices, formed by appending vertex 3 to the 
1-simplices [0, 1] and [0,2], or 



§2 = 



1 3 
2 3 



in array format. This process maybe applied recursively to develop higher dimensional simplices S3, §4, . 
as required by the application. Thus our algorithm computes simplices of the Rips complex with a hand- 
ful of sparse and dense matrix operations. 



6 Abstract Simplicial Complex Representation 

In Section 5 we saw an example of a simplicial complex which was not a manifold complex (Figure 5). 
Rips complexes described in Section 5 demonstrate one way to construct such complexes in PyDEC, 
starting from locations of vertices. There are other applications, for example in topology, where we would 
like to create a simplicial complex which is not necessarily a manifold. In addition we would like to 
do this without requiring that the location of vertices be given. For example, in topology, surfaces are 
often represented as a polygon with certain sides identified. One way to describe such an object is as an 
abstract simplicial complex [45, Section 3] . This is a collection of finite nonempty sets such that if a set is 
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in the collection, so are all the nonempty subsets of it. Figure 7 shows two examples of abstract simplicial 
complexes created in PyDEC. 




Figure 7: Examples of abstract simplicial complexes. The one of the left represents the triangulation of a Mobius 
strip and the one on the right that of a projective plane. 

In PyDEC, abstract simplicial complexes are created by specifying a list of arrays. Each array contains 

simplices of a single dimension, specified as an array of vertex numbers. Lower dimensional faces of a 
simplex need not be specified explicitly. For example the Mobius strip triangulation shown in Figure 7 
can be created by giving the array 

[0 1 3 

3 5 
3 2 5 
5 2 4 
2 4 
1 4 

as input to PyDEC. Abstract simplicial complexes need not be a triangulation of a manifold. For example 
one consisting of 2 triangles with an extra edge attached and a standalone vertex may be created using a 
list consisting of the arrays 



as input. 

The boundary matrices of a simplicial complex encode the connectivity information and can be com- 
puted from a purely combinatorial description of a simplicial complex. The locations of the vertices are 
not required. Thus the abstract simplicial complex structure is all that is required to compute these ma- 
trices as will be described in the next section. 



7 Discrete Exterior Derivative 

Given a manifold M, the exterior derivative d : Q,fiM) —>■ D.f'^^ which acts on differential p-forms, gener- 
alizes the derivative operator of calculus. When combined with metric dependent operators Hodge star, 
sharp, and flat appropriately, the vector calculus operators div, grad and curl can be generated from d. 
But d itself is a metric independent operators whose definition does not require any Riemannian metric 
on the manifold. See [1] for details. The discrete exterior derivative (which we will also denote as d) in 
PyDEC is defined as is usual in the literature, as the coboundary operator of algebraic topology [45] . Thus 

{dp a, c) = {a,dp+i c) , 
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for arbitrary p-cochain a and -chain c. Recall that the boundary operator on cochains, dp : Cp{K) —>■ 
Cp-i{K) is defined by extension of its definition on an oriented simplex. The boundary operator on a p- 
simplex = [vo,...,Vp]is given in terms of its (p - 1) -dimensional faces (p -t- 1 in number) as 



dpaf = Y,i-iy[vo,...,Vi Vp], 

i=0 



(7.1) 



where Vi means that f is omitted. Therefore, given an n-dimensional simplicial complex represented 
by §0i • • • I §w> for the discrete exterior derivative, it suffices to compute do,...,d„. As is usual in algebraic 
topology, in PyDEC we compute matrix representations of these in the elementary chain basis. Boundary 
matrices are useful in finite elements since their transposes are the coboundary operators. They are also 
useful in computational topology since homology and cohomology groups are the quotient groups of 
kernel and image of boundary matrices [45] . 

For the complex pictured in Figure 2 the boundary operators are 



do = [0 0] , 



-1 -1 

1 0-1-1 

1 0-1-1 

10 110-1 

1 1 



d2 = 






1 

-1 



1 -1 








In the following we describe an algorithm that takes as input §„ and computes both Sn-i and dn-i- This 
procedure is applied recursively to produce all faces of the complex and their boundary operator at that 
dimension. 

The first step of the algorithm converts a simplex array S„ into a canonical format. In the canonical 
format each simplex (row of §„) is replaced by the simplex with sorted indices and the relative parity of 
the original ordering to the sorted order. For instance, the simplex (1,3,2) becomes ((1,2,3),-1) since 
an odd number of transpositions (namely one) are needed to transform (1,3,2) into (1,2,3). Similarly, 
the canonical format for simplex (2, 1,4,3) is {(l,2,3,4),-l-l) since an even number of transpositions are 
required to sort the simplex indices. Since the complex dimension n is typically small (i.e. < 10), a simple 
sorting algorithm such as insertion sort is employed at this stage. We denote the aforementioned process 
caiionical_f ormat (§„) Sj^ where the rightmost column of Sjlj contains the simplex parity. Applying 
caiionical_f ormat to S2 in our example yields 
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3 




S2 = 
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3 






.2 


4 


3. 





13 1 
12 3 1 
2 3 4 -1 



Once a simplex array S„ has been transformed into canonical format, the (« - 1) -dimensional faces 
Sn-i and boundary operator d„ are readily obtained. We denote this process 



boundary_f aces(Sjlj) —>■ Sn-i,d„ . 
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In order to establish the correspondence between the n- dimensional simplices and their faces, we first 
enumerate the simplices by adding another column to S;^ to form Sj^^ . For example, 








1 


3 


1 







1 


3 


1 





s+ = 


1 


2 
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1 




1 
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3 


1 


1 




.2 


3 


4 


-1 




.2 


3 


4 


-1 


2 



The formula (7.1) is applied to S^"*" in a columnwise fashion by excluding the ?-th column of simplex 
indices, multiplying the parity column by (-1)', and carrying the last column over unchanged. For ex- 
ample. 
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1 1 
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1 2 
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1 
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1 1 


.2 
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-1 
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1 2 
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1 
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1 1 
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3 


1 2 



The resultant array is then sorted by the first n columns in lexicographical order, allowing the unique 
faces to then be extracted. 
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2 
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1 1 
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3 


-1 


1 
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3 


2 
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1 2 
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1 


1 
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3 
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1 
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3 


-1 


2 
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2 


1 1 
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1 


2 




.3 
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2 


3 


1 2 




3 


4 


-1 


2 







Furthermore, a Compressed Sparse Row (CSR) [48] sparse matrix representation of d„ as 

d„ = (ptr, indices, data) 

is obtained from the sorted matrix. For these are 

ptr=[0 1 2 3 5 7 8 9], 
indices = [0 1 1 1 2 2 2], 
data=[ 1-1 1 1-1 1-1 1 -1] 

where indices and data correspond to the fourth and third rows of the sorted matrix. 

This process is then applied to Sn-i and so on down the dimension. Since the lower dimensional 
simplex array rows are already sorted, those arrays are already in canonical format. Thus a single algo- 
rithm generates the lower dimensional faces as well as the boundary matrices at all the dimensions. The 
boundary matrices, and hence the coboundary operators, are generated in a convenient sparse matrix 
format. 
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7. 1 Generalization to Abstract Complexes 



The boundary operators and faces of an abstract simplicial complex are computed with a straightforward 

extension of the boundary_f aces algorithm. Recall from Section 6 that an abstract simplicial complex 
is specified by a list of simplex arrays of different dimensions, where the lower-dimensional simplex 
arrays represent simplices that are not a face of any higher- dimensional simplex. Generalizing the pre- 
vious scheme to the case of abstract simplicial complexes is accomplished by (1) augmenting the set of 
computed faces with the user-specified simplices and (2) modifying the computed boundary operator 
accordingly. 

Consider the abstract simplicial complex represented by the simplex arrays 



So = [5] 



§i = [l,4] 



S2 = 



1 2 
12 3 



which consists of two triangles, an edge, and an isolated vertex. Applying bouiidary_f aces to §2 pro- 
duces an array of face edges and corresponding boundary operator 



bomidary_f aces(S2) ~* Si>^2 = 



which includes all but the user-specified edge [1,4]. User- specified simplices are then incorporated into 
the simplex array in a three-stage process: (1) user-specified simplices are concatenated to the com- 
puted face array; (2) the rows of the combined simplex array are sorted lexicographically; (3) redundant 
simplices (if any) are removed from the sorted array. Upon completion, the augmented simplex array 
contains the union of the face simplices and the user-specified simplices. Continuing the example, the 
edge [1,4] is incorporated into Si as follows 
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.[1 4] 
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2 


3 



= Si 



In the final stage of the procedure, the computed boundary operator {62 in the example) is updated 
to reflect the newly incorporated simplices. Since the new simplices do not lie in the boundary of any 
higher-dimensional simplex, we may simply add empty rows into the sparse matrix representation of the 
boundary operator for each newly added simplex. Therefore, the boundary operator update procedure 
amounts to a simple remapping of row indices. In the example, the addition of the edge [1,4] into the 
fifth row of the simplex array requires the addition of an empty row into the boundary operator at the 
corresponding position. 



1 

-1 
1 








1 

-1 
1 



1 

-1 
1 










1 

-1 



1 



= d2 
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The Rips complex of Section 5 does have the location information for the vertices. However, ignoring 
those, such a complex is an abstract simplicial complex. Thus the boundary matrices for a Rips complex 
can be computed as described above. In practice some efficiency can be obtained by ordering the com- 
putation differently, so that the matrices are built as the complex is being built from the edge skeleton of 
the Rips complex. That is how it is implemented in PyDEC. 



7.2 Boundary operators and faces for cubical complexes 

The algorithm used to compute the faces and boundary operator of a given cube array (Cp Cp-i,dp) 
is closely related to the procedure discussed in Section 7 for simplex arrays. Consider a general p- 
cube in n-dimensions, denoted by the pair [{co, . . . , c„-\)ido, . . . , dp-i)] where (co,...,c„-i) are the co- 
ordinates of cube's origin and (do dp-i) are the directions which the p-cube spans. Note that the 

values [co, . . . , Cn-i,do, dp-i] correspond exactly to a row of the cube array representation introduced 
in Section 4. Using this notation, the boundary of a p-cube is given by the expression 

P-i . _ 
dp[{co,...,Cn-i]{do,...,dp-i)] = (-l)'([(co,...,Cd. + 1 Cn-i]{do,...,di,...,dp-i)]- 



i=0 



[(Co, ... , -I- 0, ... , Cn-i){do, ...,di,..., dp-i)]) (7.2) 



where di denotes the omission of the i-th spanning direction and q,. is the corresponding coordinate. 
For example, the boundary of a square centered at the location (10, 20) is 



a2[(10,20)(0,l)] = [(11,20)(1)] - [(10,20)(1)] - [(10,21)(0)] + [(10,20)(0)]. 



(7.3) 



The canonical format for a p-cube is the one where the spanning directions are specified in ascend- 
ingorder. For instance, the 2-cube [(10,20)(0, 1)] is in the canonical format because (io < c?i . As with sim- 
plices, each cube has a unique canonical format, through which duplicates are easily identified. Since 
the top-level cube array C„ is generated from a bitmap it is already in the canonical format and no re- 
ordering of indices or parity tracking is necessary. 

Applying Equation 7.2 to a p-cube array with N members generates a collection 2N oriented faces. 
In the mesh illustrated in Figure 3 the three squares in Cz are initially expanded into 



C2 = 





1 

1 1 











1 

1 1 

1 

2 
2 1 



1 
1 1 

1 

1 1 
1 2 












where the fourth column of encodes the orientation of the face and the fifth column records the 2- 
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cube to which each face belongs. Sorting the rows of in lexicographical order 
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1 2 
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1 1 


2 









= Ci 



allows the unique faces to be extracted. Lastly, a sparse matrix representation of the boundary operator 
is obtained from the sorted cube array in the same manner as for simplices. 



8 Review of Whitney Map and Whitney Forms 

In this section we review and collect some material, most of which is well-known in DEC and finite ele- 
ment exterior calculus research communities. It is included here partly to fix notation. In this section we 
also give the monomials based definition of inner product of differential forms. This is not the way inner 
product of forms is usually defined in most textbooks, [44] being one exception we know of. The mono- 
mial form leads to an efficient algorithm for computation of stiffness and mass matrices for Whitney 
forms given in Section 9. 

The basic function spaces that are useful with exterior calculus are the space of square integrable p- 
forms on a manifold and Sobolev spaces derived from that. Let M be a Riemannian manifold, a manifold 
on which a smoothly varying inner product is defined on the tangent space at each point. Let g be its 
metric, a smooth tensor field that defines the inner product on the tangent space at each point on M. 

For differential forms on such a manifold M, the space of square integrable forms is denoted L^Q.P (M) . 
One can then define the spaces HDP{M) which generalize the spaces H{Aiv) and //(curl) used in mixed 
finite element methods [4] . To define I?D}' {M) one has to define an inner product on the space of forms 
which is our starting point for this section. All these function spaces have been discussed in [2] and [4]. 
The definitions and properties of Whitney map and Whitney forms is in [21], the geometric analysis 
background is in [38] and the basic definition of inner products on forms is in [1, page 411]. 

8. 1 Inner product of forms 

To define the spaces L^D-P^M) and HDP{M) more precisely, we recall the definitions related to inner 
products of forms. We will need the exterior calculus operators wedge product and Hodge star which 
we recall first. For a manifold M the wedge product a : £if{M) x f2^(M) — ► D.P^'^iM) is an operator for 
building ip + q) -forms from p-forms and (/-forms. It is defined as the skew-symmetrization of the tensor 
product of the two forms involved. For a Riemannian manifold of dimension n, the Hodge star operator 
* : QPiM) — Q."~P{M) is an isomorphism between the spaces of p and {n - p)-forms. For more details, 
see [1, page 394] for wedge products and [1, page 411] for Hodge star. 
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Definition 8.1. Given two smooth p-forms a, )6 £ QP(Af) on a Riemannian manifold M, their pointwise 
inner product at point x £ M is defined by 

(a;(jc),j6(x)>ju = a(x)A*j6(x), (8.1) 

where jU = * 1 is the volume form associated with the metric induced by the inner product on M. 

The pointwise inner products of forms can be defined in another way, which will be more useful to 
us in computations. The second definition given below in Definition 8.2 is equivalent to the one given 
above in Definition 8.1. The operator (j (the sharp operator) used below is an isomorphism between 1- 
forms and vector fields and is defined by g(a;'*,X) = a (X) for given 1 -form a and all vector fields X. See [1] 
for details. 

Definition 8.2. Let ai,...,ap and j6i, . . . , )3p be 1-forms on a Riemannian manifold M. By analogy with 
polynomials we'll call p-forms of the type aiA--- Aap and )6i A • • • A )6p monomial p -forms. Define the 
following operator at a point x £ M: 

<ai A • • • A ap, )6i A • • • A )8p) := det[g(a;., iSj.) ] , (8.2) 

where is the matrix obtained by taking 1 < i,j < p. In the equation above, all the 1-forms 

are evaluated at the point x £ M. Extend the operation in (8.2) bilinearly pointwise to the space of all 
p-forms. It can be shown that this defines a pointwise inner product of p-forms equivalent to the one 
defined in (8.1). Note that if a, = Pi for all i, the expression on the right in (8.2) is the Gram determinant. 

Remark 8.3. Note that unlike (8.1) the definition in (8.2) does not involve wedge product and Hodge 
star explicitly. This is an advantage of the latter form since a discrete wedge product is not available in 
PyDEC. The RHS of (8.2) does involve the sharp operator, but as we will see in the next section, this is 
easy to interpret for the purpose of discretization in this context. 

Definition 8.4. The pointwise innerproduct in (8. 1), or equivalentiy in (8.2), induces an inner product 
on Mas 

ia,Ph^= [ {aM,Pix))p. (8.3) 
Jm 

The space of p-forms obtained by completion of QP(M) under this inner product is the Hilbert space of 
square integrable p-forms L^D.P{M). The other useful space mentioned at the beginning of this section 
is the Sobolev space HOP^M) := {a £ L^D.P{M) \ da £ L^nP+^{M)]. 

8.2 Whitney map and Whitney forms 

Let K be an n-dimensional manifold simplicial complex embedded in IR^ and \K\ the underlying space. 
The metric on the interiors of the top dimensional simplices of K will be the one induced from the em- 
bedding Euclidean space R^. As is usual in finite element methods, finite dimensional subspaces of the 
function spaces described in the previous paragraph are used in the numerical solution of PDEs. The 
finite dimensional spaces can be obtained by "embedding" the space of cochains into these spaces by 
using an interpolation. For example, to embed C^(iC;IR) into L^D.P[\K\] one can use the Whitney map 
W : CP{K;U) L^QP{\K\), which will be reviewed in this subsection. The image W(CP(K;K)) is a linear 
vector subspace of L^D.p[\K\) and is the space of Whitney p-forms [13, 21, 52] and is denoted ^[D.P{\K\) 
in [3, 4]. (We use instead of used in [4].) The embedding of cochains is analogous to how scalar val- 
ues at discrete sample points would be interpolated to get a piecewise affine function. In the scalar case 
also, the space of such functions is a vector subspace of square integrable functions. In fact, W(C^(^r; R)), 
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the space of Whitney 0-forms is the space of continuous piecewise affine functions on\K\. The Whitney 
map for p > is actually built from barycentric coordinates which are the building blocks of piecewise 
linear interpolation. Thus the embedding of Cf{K;U) into L^QP(|Ar|) for p > can be considered to be 
a generalization of the embedding of C" {K; U) into L^DP[\K\) . Thus Whitney forms enable only low or- 
der methods. However, arbitrary degree polynomial spaces suitable for use in finite element exterior 
calculus have been discovered [2, 3]. These however are not yet a part of PyDEC. 

The space of Whitney p-forms is the space of piecewise smooth differential p-forms obtained by ap- 
plying the Whitney map to p-cochains. It can be thought of as a method for interpolating values given 
on p-simplices of a simplicial complex. For example, inside a tetrahedron Whitney forms allow the inter- 
polation of numbers on edges or faces to a smooth 1-form or 2-form respectively. As mentioned above, 
for 0-cochains, i.e. scalar functions sampled at vertices, the interpolation is the one obtained using the 
standard scalar piecewise affine basis functions on each simplex, that is the barycentric coordinates cor- 
responding to each vertex of the simplex. We recall the definition of barycentric coordinates followed by 
the definition of the Whitney map. 

Definition 8.5. Let af = [vq,..., Vp] be a p-simplex embedded in K^. The affine functions p, : IR^ K, 
i = 0,...,p, which when restricted to cr^ take the value 1 on vertex Vi and on the others, are called the 

barycentric coordinates in cr^. 

Definition 8.6. Let be an oriented p-simplex [f ((,, . . . , in an ^-dimensional manifold complex K, 
and [a^Y the corresponding elementary p-cochain. We define 

P 

W((o-P)*) := p! X (-1) Mifc dfii^ A • • • A dp,-, a • • • a djui^ , (8.4) 

fc=0 

where p/j. is the barycentric coordinate function with respect to vertex y,-, and the notation dp,-, indicates 
that the term rfp;, is omitted from the wedge product. The Whitney map W: CP{K;U) L^D.P[\K\) is 
the above map W extended to all of CP{K;U) by requiring that W be a linear map. W{(_aP)*) is called 
the Whitney form corresponding to a^, and for a general cochain c, W(c) is called the Whitney form 
corresponding to c. 

For example, the Whitney form corresponding to the edge [ f o , f i ] is W( [ , f i ] * ) = po d pi - pi d po , 
and the Whitney form corresponding to the triangle [i^i, 1^2, i^s] in a tetrahedron [I'o. 1^1. i'2> 1*3] is 

W([i;i, V2, V3]*) = 2 (pi dp2 A dp3 - p2 dpi A dps + Pa dpi a dp2) . 

Remark 8.7. If we were using local coordinate charts on a manifold then at any point a p-form would be 
a linear combination of monomials. Note from (8.4) that the Whitney form W(c7*) is a sum of monomi- 
als with coefficients. Thus Whitney forms allow us to treat forms at a point as a linear combination of 
monomials even though we are not using local coordinate charts. 

We emphasize again that this section was a review of known material. We have tiled to present this 
material in a manner which makes it easier to explain the examples of Section 1 1 and the construction 
of mass matrix for Whitney forms described in the next section. 

9 Whitney Inner Product of Cochains 

Given a manifold simplicial complex K, an inner product between two p-cochains a and b can be de- 
fined by first embedding these cochains into L^D,f[\K\) using Whitney map and then taking the inner 
product of the resulting Whitney forms [21] . 
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Definition 9.1. Given two p-cochains a, be Cf{K;U), their Whitney inner product is defined by 



J\K\ 

using the inner product on forms given in (8.3). The matrix for Whitney inner product of jo-forms 
in the elementary p-cochain basis will be denoted Mp. That is, Mp is a square matrix of order Np (the 



number i and j respectively. 

The integral in (9.1) is the sum of integrals over each top dimensional simplex in K. Inside each such 
simplex the inner product of smooth forms applies since the Whitney form in each simplex is smooth 
all the way up to the boundary of the simplex. The interior of each top dimensional simplex is given an 
inner product that is induced from the standard inner product of the embedding space K^. 

Remark 9.2. Given cochains a,beC'^{K;M)we will refer to their representations in the elementary cochain 
basis also as a and b. Then the matrix representation of the Whitney inner product of a and b is (a, b) = 
a^Mpb. 

The inner product of cochains defined in this way is a key concept that connects exterior calculus 
to finite element methods and different choices of the inner product lead to different discretizations 
of exterior calculus. This is because the inner product matrix Mp is the mass matrix of finite element 
methods based on Whitney forms. The details of the efficient computation of Mp for any p and n will be 
given in Section 9.2 and 9.3. 

Recall that for a Riemannian manifold M, if 5p+i : nP+i(M) — nP{M) is the codifferential, then the 
Laplace-deRham operator on p-formsis Ap := dp-i5p + 5p+i dp. For a boundaryless M the codifferential 
6p+i is the adjoint of the exterior derivative dp. In case M has a boundary, we have instead that 



See [1, Exercise 7.5E] for a derivation of the above. Now consider Poisson's equation ApU = fon p-forms 
defined on a p-dimensional simplicial manifold complex K. For simplicity, we'll consider the weak form 
of this using smooth forms rather than Sobolev spaces of forms. See [2, 4] for a proper functional analytic 
treatment. We will also assume that the correct boundary conditions are satisfied, so that the boundary 
term in (9.2) is 0. In our simple treatment, the weak form of the Poisson's equation is to find awe (|i<r|) 
such that (dp u, 5p v)i2 + (dp u, dp v)i2 = (/, v)i2 . Thus it is clear that a Galerkin formulation using Whit- 
ney forms ^j"QP will require the computation of a term like (dp Wa, dpWb)i2 for cochains a and b. By 
the commuting property of Whitney forms dpW = W dp (where the second dp is the coboundary opera- 
tor on cochains) we have that the above inner product is equal to (W dp a, Wdp b)i2. (See [21] for a proof 
of the commuting property.) Now, by definition of the Whitney inner product of cochains in (9.1) this 
is equal to (dp a, dp b) in the inner product on (p -I- 1) -cochains. The matrix form of this inner product 
can be obtained from the mass matrix Mp+i as d^ Mp+i dp. This is what we mean when we say that the 
stiffness matrix can be computed easily from the mass matrix. The term on the right in the weak form 
wiU use the mass matrix Mp. Since codifferential of Whitney forms is 0, the first term in the weak form 
has to be handled in another way, as described in [4] . 

Remark 9.3. Exploiting the aforementioned commutativity of Whitney forms to compute the stiffness 

matrix represents a significant simplification to our software implementation. While computing the 
stiffness matrix directiy from the definition is possible, it is a complex operation which requires consid- 
erable programmer effort, especially if the performance of the implementation is important. In contrast. 




(9.1) 





(9.2) 
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our formulation requires no additional effort and has the performance of the underlying sparse matrix- 
matrix multiplication implementation, an optimized and natively- compiled routine. All of the complex 
indexing, considerations of relative orientation, mappings between faces and indices, etc. is reduced to 
a simple linear algebra expression. The lower dimensional faces are oriented lexicographically and the 
orientation information required in stiffness matrix assembly is implicit in the boundary matrices. 

9. 1 Computing barycentric differentials 

Given that the Whitney form W(c7*) in (8.4) is built using wedges of differentials of barycentric coordi- 
nates, it is clear that the algorithm for computing an inner product of Whitney forms involves compu- 
tation of the gradients or differentials of the barycentric coordinates. The following lemma shows how 
these are computed using simple linear algebra operations. 

Lemma 9.4. Let = [vo,...,Vp] be a p -simplex embedded in , p < N where the vertices Vi £ are 
given in some basis for . Let XeU^^P bea matrix whose j -th column consists of the components ofd p,j 
in the dual basis, for j = l,...,p . Let Vq e U^"" ^ be a matrix whose j-th column is vj - vo, for j = l,...,p. 
Then = {VqVo)~ Vq = , the pseudoinverse ofV^. 

Proof. Let C = [/Ji , . • . , /ip] ^ be the vector of barycentric coordinates (other than /Jq) with respect to a^, for 

a point X = [xi x^] ^ e IR^. Then by definition of barycentric coordinates and simplices, Vo( = x - vq 

is the linear least squares system for the barycentric coordinates. Thus ( = [x-vo) which implies that 
d( = = V+. □ 

Remarks. 5. The use of normal equations in the solution of the least squares problem in the above propo- 
sition suffers from the well-known condition squaring problem. This is only likely to be a problem if the 

simplices are nearly degenerate. In that case one can just use an orthogonalization method to compute a 
QR factorization and use that to solve the least squares problem. Notice that in typical physical problems 
Vq will typically be 2 x 2, 3 x 2 or 3 x 3 matrix so any of these methods are easy to implement. 

Once the components for dpi have been obtained for i = l,...,p, the components of dpo can be 
obtained by noting that djuo -i- • • • -i- d jUp = which follows from the fact that the barycentric coordinates 
sum to 1. Also note that the components of the gradients Vpi will be the same as those of dpi if the 
standard metric of Euclidean space is used for the embedding space which is the case in all of PyDEC. 

9.2 Whitney inner product matrix 

We will now use the inner product of forms in (8.2) and the cochains inner product defined in (9.1) to 
give a formula for the computation of Mp, the Whitney inner product matrix for p-cochains. We will also 
refer to this as the Whitney mass matrix. 

Notation 9.6. Given simplices a and t the notation a >t ovt <a means t is a face of a. Note that this 
means t can be equal to a since any simplex is its own face. For proper inclusion we use t < a or a> r 
to indicate that t is a proper face of a. The use of this notation simplifies the expression of summations 
over various classes of simplices in a complex. For example, given two fixed p-simplices af and 

L 

a" 

' J 

is read as "sum over all ra-simplices a" which have and as faces" . Another notation used in the 
proof of the proposition below is the star of a simplex a, written St(a) (not to be confused with the dual 
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star ★fj). This star St((r) is the union of the interiors of all simplices of the complex that have cr as a face. 
That includes a also. The closure of this open set is called the closed star and written Sta. This is the 
union of simplices that contain a. 

Proposition 9.7. Leta^ = [Vig vt^^] anda^ = [vj^,. . . , vj^] be oriented p-simplices in an n-dimensional 

manifold simplicial complex K, with < p < n. Then the row i, column j entry of the Whitney inner 
product matrix Mp is given by 



Mp{i,j] = {p^f £ Zi-lf^'ckJ lii,lij,p, 
a" k,l=Q •'o'" 



' J 



where Cki = 1 for p = 0, and for p>0 

<d)U,„,d)U;„> ... {dpi^,dpji) ... <d)U,„,d)U;p> 



Cki = det 



id Pi,, dp j,} 



{dpi,,dpj,} 



{dpi^,dpj,} ... {dpi^.dpj^} ... {dpi^,dpj^) 

the hats indicating the deleted terms. Here p is the volume form corresponding to the standard inner 
product in and pt, and pj, are the barycentric coordinates corresponding to vertices ik and ji . 

Proof See Appendix. □ 

For the p = n case a simpler formulation is given in Proposition 9.9. The above proposition shows 
that computation of the Whitney inner product matrix involves computations of inner products of dif- 
ferentials of barycentric coordinates. Since the only metric implemented in PyDEC is the standard one 
inherited from the embedding space K^, 

{dpi,dpj) = g{{dpi)K {dpj)h = Vpi ■ Vpj . 

Example 9.8. Consider the simplicial complex corresponding to a tetrahedron embedded in for 
which we want to compute M2, the Whitney inner product matrix for 2-cochains. Here N = n = 3 and 
p = 2 and M2 is of order N2 = 4, the number of triangles in the complex. Label the vertices as 0, 1,2,3. 
Then by PyDEC's lexicographic numbering scheme, the edges numbered to 5 are [0,1], [0,2], [0,3], [1,2], 
[1, 3] , [2, 3] and the triangles numbered to 3 are [0, 1,2], [0, 1, 3] , [0, 2, 3], and [1, 2, 3] . We will describe the 
computation of the row 0, column 3 entry and the row 0, column 1 entry of M2. The (0,3) entry corre- 
sponds to the inner product of cochains [vo,Vi,V2]* and [1^1,^2,^3]* since, in the lexicographic ordering 
and naming convention of PyDEC these are and (T3 respectively. Thus we are computing 

The corresponding Whitney forms are 

W(o-o)* = 2! [po dpiAdp2- pidpoA dp2 + p2dpoA dpi) 
MV[al)* =2\{pi dp2Adp3- p2dpiAdp3 + p3dpiAdp2). 

Then {W{aiy,W{aiy)^J{2\f is 

/ jUojUi(djUi AdjU2,d)U2 AdjU3))U- / jUoAi2(djUi Ad)U2,djUi Ad)U3)jU-l-..., (9.3) 
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where ju is just dx a dy a dz, the standard volume form in K^. Each term like (djUi AdjU2,d)U2 AdjUs) is 



det 



<djUi,d)U2> 

(djU2,d)U2) 



<djUi,d)U3> 

(djU2,d)U3) 



which in the notation of Prop. 9.7 is 



det 



<djZi^d/Ji> 

<d/Z2,djLli> 



<d;U0,djU2> 

<d;Ui,d^2> 

<d)U2,djU2> 



<d;U0,djLi3> 
<d;Ul,d^3> 
(d|U2,djLi3> 



Using a shorthand notation for matrices like above, the 2x2 matrices whose determinants need to be 
computed for calculating the (0, 3) entry of M2 are given below. 
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11 12 13 
01 02 03 
21 22 23 



12 11 13 
02 01 03 
22 21 23 



13 11 12 
03 01 02 
23 21 22 



21 22 23 
01 02 03 
n 12 13 



22 21 23 
02 01 03 
12 11 13 



23 21 22 
03 01 02 
13 11 12 



Removing the deleted rows and columns the above matrices are given below as the actual 2x2 matrices. 
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(9.4) 



The 2x 2 matrices whose determinants are needed in computing {{(^D* ,{(^l]*]> i-e., entry (0,1) of M2 are 
given below. 
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(9.5) 



Recall that each number in these matrices is shorthand for an inner product of two barycentric dif- 
ferentials. For example, the entry 12 stands for (djUi, d;U2> = g((d)Ui)^ (d;U2)') = VjUi • VjU2- 
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Proposition 9.9. For p = n, Mp is a diagonal matrix with Mp{i, i) = l/|a"|, where \a"\ is the volume of 
the simplex. 

Proof. For any «-simplex it", the Whitney form W{cr") * is on other n-simplices and so Mp is diagonal. 
Furthermore, it is a constant coefficient volume form on a" with f^„W{a")* = 1. See [21] for proofs of 
these properties. Thus it must be that W(c7") * = i^/\a"\ where ju is the volume form on the simplex. Thus 

ma") * , W(a") *)n = W{a") * a * W((t") * = . 
Thus (W((t")*,W((t")*)^2 is/^„ju/|o-"l^ whichis l/\a"\. □ 



9.3 Algorithm for Whitney inner product matrix 

We motivate our algorithm for Whitney mass matrix computation by making some observations about 
Example 9.8. The first, and obvious observation is that matrix Mp is symmetric, being an inner product 
matrix. Thus only the diagonal entries and those above (or below) the diagonal need be computed. A 
more interesting efficiency comes from the structure of the entries of the matrix collections, such as ones 
shown in (9.4) and (9.5). Note that many entries repeat in the shorthand collection of matrices in (9.4) 
and (9.5). For example the entry 12 appears 4 times by itself in the matrix collection (9.4). Moreover, 
due to the symmetry of inner product, the entry 12 corresponds to the same result as the entry 21 and 
21 appears 4 times as well. That entry also appears 4 times in the collection (9.5). Thus it is clear that a 
saving in computational time can be achieved by doing such calculations only once. That is, <d jUi , d jU2> = 
(d)U2, djUi) need only be computed once for the tetrahedron. 

The determinants of all the matrices in a collection such as (9.4) are needed to plug into an expression 
like (9.3) to obtain a single entry (in this case row 0, column 3) of the Whitney inner product matrix for 
p-cochains (p = 2 in this case), whose size (4 x 4 in this case) depends on the number of p-simplices in 
the simplicial complex. Thus reusing repeated inner products of barycentric differentials can add up to 
a substantial saving in computational expense when all the unique entries of Mp are computed. These 
savings are quantified later in this subsection. 

Another useful point to note in the example calculation is that the collection (9.5) of matrices can be 
obtained from the collection (9.4) by keeping the first digit in each entry same and making the substi- 
tutions 1 ^ 0; 2 ^ 1; and 3 — ' 3 in the second digit. The first digits in the two collections are the same 
because both correspond to the triangle cTq. The substitution above works for the second digit because 
cjg = [vi,V2,V3] and af = [vo,vi,V3]. This suggests the use of a template simplex for creating a template 
collection of matrices whose determinants are needed. The actual instances of the collections can then 
be obtained by using the vertex numbers in a given simplex. This is another idea that is used in the 
algorithm implemented in PyDEC. The algorithm takes as input a manifold simplicial n-complex K, em- 
bedded in and < p < n. The output is Mp, an Np x Np matrix representation of inner product on 
CiK; M) using elementary cochain basis. If a naive algorithm, which does not take into account the du- 
plications in determinant calculations were to be used, the number of operations required in the mass 
matrix calculation are 



iV„x 







'n + l\ 









n 



2 



xNp^xiOipl) orO(p^)). 



The last term is written as 0{p\) or 0[p^) because a determinant can be computed using the formula for 
determinant or by LU factorization. For low values of p (i.e. < about 5) the formula will likely be better. 
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According to the above formula, for example, for n = 3,p = 2, the number of determinants required 
in a naive implementation of mass matrix calculation would be 



+ 



= 10x9 = 90. 



But there are only 21 unique determinants needed for n = 3,p = 2. Our algorithm computes the unique 
determinants first and the operation count is 



iV„x 



(n+l 
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'n + l 
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[ P , 
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X Np^ X (0(p!) or 0(p^)) . 



Figure 8 shows a comparison of determinant counts for our algorithm compared with a naive al- 
gorithm that does the duplicate work that our PyDEC algorithm avoids. Note that for any n, the most 
advantage is gained intermediate values of p. The savings that the PyDEC implementation provides over 
a naive algorithm are several orders of magnitude, especially for moderately large n and higher. For p = n 
case, in PyDEC we use the shortcut described in Proposition 9.9. 



10 Metric Dependent Operators 

We now describe the PyDEC implementations of some metric dependent exterior calculus operators. 

The simplicial complex K is now supposed to be an approximation of a Riemannian n-manifold M. 
The metric implemented in PyDEC is the one induced from an embedding space IR^. The main metric 
dependent operator is the Hodge star which enables the discretization of codifferential and Laplace- 
deRham operators. The sharp and flat, which are isomorphisms between 1 -forms and vector fields, are 
not implemented. 

For the DEC Hodge star, the implementation is using the circumcentric dual as in [19, 31] and the 
other operators are then simply defined in terms of the exterior derivative and the Hodge star. For Py- 
DEC's implementation of low order finite element exterior calculus, we define the Hodge star to be the 
Whitney mass matrix described in Section 9. The other operators are defined by analogy with DEC even 
though the dual mesh concept is not part of finite element exterior calculus. Extensive experimental 
justification for this approach can be seen in its effectiveness in numerical experiments in [8] and in [32] . 

For the definitions in this section we will need two cochain complexes of real-valued cochains. One 
will be on the simplicial complex K and for brevity we'll call this space of p-cochains CPiK) instead of 
C^{K;U). The other cochain complex is on the circumcentric dual cell complex -kK and we'll denote 
the in - p) -dimensional cochains as D"~P[-k K). At each dimension, these will be connected by discrete 
Hodge star operators to be defined below. Since the exterior derivative is the coboundary operator, the 
matrix representation for the exterior derivative on the dual mesh is the boundary operator. The matrix 
form for the DEC Hodge star on p-cochains will be denoted *p : CiK) — ► D'^~P{-kK). One box of the 
primal and dual complexes is shown below. 
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Figure 8: Comparison between PyDEC and a naive algorithm for computing Whitney mass matrix. The figure 
shows the number of determinant computations needed by the two algorithms, for various values of p, the Whit- 
ney form dimension, and n, the simplicial complex dimension. The embedding dimension is not relevant in these 
calculations. For p-n case we use the shortcut described in Proposition 9.9 so that case is not shown. 
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As described in [19, 31] and other references, the DEC Hodge star is defined by 

{*pa*,-kaj) ^ {a*,aj) 
\*aj\ \aj\ 

for p-simplices a; and cry. Here a* is the elementary cochain corresponding to a and {a* ,r) stands 
for the evaluation of the cochain a* on the elementary chain t. Thus the matrix representation of the 
DEC Hodge star *p is as a diagonal matrix with *p(j, i) = \-kOi\l\Oi\. In [19, 31] this was defined for well- 
centered meshes. For the codimension 1 Hodge star the definition extends to Delaunay meshes with a 
slight additional condition for boundary simplices. This extension involves computing the volume of 
•k a taking into account signs. Consider a codimension 1 simplex a shared by simplices L and R. For 
the portion of ★ cr corresponding to L, the sign is positive if the circumcenter and remaining vertex of 
L are on the same side of a. Similarly for R. (For surface meshes and higher dimensional analogs the 
circumcenter condition above is one way to define a Delaunay-like condition.) If cr is a codimension 1 
face of top dimensional t and is on domain boundary then the circumcenter of t and vertex opposite 
to a should be on the same side. The smooth Hodge star on p-forms satisfies * * = {-1)P^"'~p\ In the 
discrete setting * * is written as *~i *p or *p*p^ and this is defined to be (-l)P("~P' / where I is the 
identity matrix. 

In the smooth theory, the codifferential <5p+i : QP+i(Af) - nP[M) is defined as dp+i = (-1)"P+i * d* 
and so we define the discrete codifferential Sp+i : CPiK) - C"-P(i<0 as 6p+i := (-l)"P+i *-id^*p+i. 
For finite element exterior calculus implemented in PyDEC, we take this to be the definition, without 
reference to a dual mesh. If we now take * p to be the Whitney mass matrix then dp and 5p+ 1 are adjoints 
(up to sign) with respect to the Whitney inner product on cochains as shown in [32] . We will call the use 
of Whitney mass matrix as * p to be a Whitney Hodge star mattix. 

In the discrete setting the Laplace- deRham operator is implemented in the weak form. For 0< p< n 
the discrete definition is Ap := dp *p+i dp + (-1)'''"^''""''"'"^' *p dp-i dp_j *p, with the appropriate 
term dropped for the p = and p = n cases. The above expression involves inverses of the Hodge star, 
which is easy to compute for DEC Hodge star since that is a diagonal matrix. For a Whitney Hodge 
star see [8, 32] for various approaches to avoiding explicitiy forming the inverse Whitney mass matrix in 
computations. 



10.1 Circumcenter Calculation 

Circumcentric duality is used in DEC. To compute the DEC Hodge star, a basic computational step is the 
computation of the circumcenter of a simplex. We give here a linear system for computing the circum- 
center using barycentric coordinates. 

The circumcenter of a simplex is the unique point that is equidistant from all vertices of that simplex. 
In the case that a simplex (or face) is not of the same dimension as the embedding (e.g. a triangle em- 
bedded in M*) , we choose the point that lies in the affine space spanned by the vertices of the simplex. In 
either case we can write the circumcenter in terms of barycentric coordinates of the simplex. 

Let (jP be the p-simplex defined by the points {vo,Vi,... Vp} in U^. Let R denote the circumradius 
and c the circumcenter of simplex, which can be written in barycentric coordinates as c = bj vj where 
bj is the barycentric coordinate for the circumcenter corresponding to vj. For each vertex / we have 

2 



P 
7=0 



-i?^ = 0. 



which can be rewritten as 



Vi-Vi-2vi 



P 

\j=o 



p 

7=0 



■R^ = 0. 
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Here the norm and the dot product are the standard ones on U". Rearranging the above yields 



\J=o ) 



p 



The second term on the left hand side is some scalar which is unknown, but is the same for every equa- 
tion. So we can replace it by the unknown Q and write 



2vi 



Y.bjVj\ + Q=Vi-Vi. 



With the additional constraint that barycentric coordinates sum to one, we have a linear system with 
p + 2 unknowns (&o • • • bp and Q) and p + 2 equations with the following matrix form 



2vi-vo 2vi-vi 



2vp-vo 2vp-v\ 
1 1 
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2vi-Vp 1 
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The solution to this yields the barycentric coordinates from which the circumcenter c can be located. 
Another quantity required for DEC Hodge star is the unsigned volume of a simplex. This can be com- 
puted by the well-known formula VdetV^V I pi where V is the phy N matrix with rows formed by the 
vectors {i^i - vo, V2-vo,...Vp- vq}. 



1 1 Examples 

In the domains for which PyDEC is intended, it is often possible to easily translate the mathematical 
formulation of a problem into a working program. To make this point, and to demonstrate a variety of 
applications of PyDEC, we give 5 examples from different fields. The first example (Section 11.1) is a res- 
onant cavity eigenvalue problem in which Whitney forms work nicely while the nodal piecewise linear 
Lagrange vector finite element fails when directly applied. The second is Darcy flow (Section 11.2), 
which is an idealization of the steady flow of a fluid in a porous medium. We solve it here using DEC. The 
third problem (Section 11.3) is computation of a basis for the cohomology group of a mesh with several 
holes. This is achieved in our code here by Hodge decomposition of cochains, again using DEC. Next 
example (Section 11.4) is an idealization of the sensor network coverage problem. Some randomly lo- 
cated idealized sensors in the plane are connected into a Rips complex based on their mutual distances. 
Then a harmonic cochain computation reveals the possibility of holes in coverage. The last example 
(Section 1 1.5) involves the ranking of alternatives by a least squares computation on a graph. 

None of these problem is original and they all have been treated in the literature by a variety of tech- 
niques. We emphasize that we are including these just to demonstrate the capabilities of PyDEC. We have 
included the relevant parts of the Python code in this paper. The full working programs are available with 
the PyDEC package [7]. 



11.1 Resonant cavity curl-curl problem 

An electromagnetic resonant cavity is an idealized box made of a perfect conductor and containing no 
enclosed charges in which Maxwell's equations reduce to an eigenvalue problem. Several authors have 
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popularized this example as one of many striking examples that motivate finite element exterior cal- 
culus. See for instance [4]. The use of S^^ finite element space, i.e. piecewise linear, Lagrange finite 
elements with 2 components, yields a corrupted spectrum. On the other hand, the use of elements, 
i.e., Whitney 1-forms yields the qualitatively correct spectrum. For detailed analysis and background 
see [4, 12]. 

Let M c be a square domain with side length n. We first give the equation in vector calculus 
notation and then in the corresponding exterior calculus notation. In the former, the resonant cavity 
problem is to find vector fields E and eigenvalues A £ K such that 

curl curl £■ = onM and £'||=0on5M, 

where is the tangential component of E on the boundary. Here cnrXcp = {dcp/dy, -dcpldx) and curl v = 
dv2ldx-dv\ldy for scalar function (p and vector field v = (i^i, 1^2)- Note that for A ^ this equation is 
equivalent to the pair of equations AE = XE and divE = 0. This is because the vector Laplacian A = 
curl o curl - grad o div and divo curl = 0. 

Now we give the equation in exterior calculus notation so the transition to PyDEC will be easier. Let 
u £ D}{M) be the unknown electric field l-form and i : dM •—>■ M the inclusion map. Then the above 
vector calculus equation is equivalent to 

§2 di M = A M in M 
i*u = ondM. 

The puUback /*« by inclusion map means restriction of u to the boundary, i.e., allowing only vectors 

tangential to the boundary as arguments to u. As usual, we will seek u not in f2^(M) but in HD.^{M) 
subject to boundary conditions. Define the vector space V = {v\ve HD.^(_M), i* v = Oon dM}. 

To express the PDE in weak form, we seek a (a. A) in x R such that {^2 di u, v)i2 = X{u, v)i2 for all £ 
V. By the properties of the codifferential, the expression on the left is equal to (di u, di v)i2 -fgj^ ua * di v. 
But the boundary term is because u is in V. Thus the weak form is to find a (u. A) £ V x R such that 
(di M,di v)i2 = A(m, v)i2 for all veV. 

Taking the Galerkin approach of looking for a solution in a finite dimensional subspace of V here we 
pick the space of Whitney 1-forms, that is, ^^D,^ as the finite dimensional subspace. We define these 
over a triangulation of M which we will call K. The Whitney map W : C^(K;[R) L^n^[\K\) is an injection 
with its image ^^{K). Thus an equivalent formulation is over cochains. Using the same names for the 
variables, we seek a (m,A) £ C^iK;U) x R such that (di W m, di W i;)i2 = A(W u,Wv)i2 for all 1-cochains 
V £ C^(K;R). Since the Whitney map commutes with the exterior derivative and coboundary operator, 
and using the definition of cochain inner product, the above is same as (di u, di v) = X{u, v) where now 
the inner product is over cochains and di is the coboundary operator. In matrix notation, using *i and 
*2 to stand for the Whitney mass matrices Mi and M2, the generalized eigenvalue problem is to find 
(m. A) £ {K; R) X R such that 

df *2 di M = A *i M. 

We now translate this equation into PyDEC code. Once the appropriate modules have been im- 
ported, a simplicial complex object sc is created after reading in the mesh files. Now the main task is to 
find matrix representations for the stiffness matrix df *2 di and the mass matrix * 1 . This is accomplished 
by the following two lines, where K is the stiffness matrix : 

K = sc[l].d.T * whitney_innerproduct (sc , 2) * sc[l].d 
M = wliitney_innerproduct (sc , 1) 

The boundary conditions can be imposed by simply removing the edges that lie on the boundary. The 
indices of such edges is easily determined and stored in the list iioii_bo\mdary_iiidices which is used 
below to impose the boundary conditions : 
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Figure 9: The first 50 nonzero eigenvalues for the resonant cavity problem of Section 11.1 and the eigenvector 
corresponding to one of these eigenvalues. The eigenvector is a 1-cochain which is visualized as a vector field by 
first interpolating it using a Whitney map. See Section 11.1 for details. 

K = K [non_boundary_indices ,:][:, non_boundary_indices] 
M = M [iioii_boundary_indices ,:][:, non_boundary_indices] 

Now all that remains is to solve the eigenvalue problem. To simplify the code and because the matrix size 
is small, we use the dense eigenvalue solver scipy . linalg . eig 

eigenvalues, eigenvectors = eig (K . todense () , M.todenseO) 

Some of the resulting eigenvalues are displayed in the left part of Figure 9. The 1-cochain u which is 
the eigenvector corresponding to one of these eigenvalues is shown as a vector field in the right part 
of Figure 9. The visualization as a vector field is achieved by interpolating the 1-cochain u using the 
Whitney map and then sampling the vector field (W w)" at the barycenter. This is achieved by the PyDEC 
command: 

bases, arrows = s iinplex_quivers ( sc , all_values ) 

where all_values contains both the known and the computed values of the 1-cochain. There is no 
sharp operator in PyDEC. But since PyDEC only implements the Riemannian metric from the embedding 
space of simplices the transformation from 1-form to vector field just involves using the components of 
the Whitney 1-form as the vector field components. 

11.2 Darcy flow or Poisson's in mixed form 

We give here a brief description of the equations of Darcy flow and their PyDEC implementation. For 
more details see [34] . The resonant cavity example in Section 11.1 was implemented using finite element 
exterior calculus. For variety we use a DEC implementation for Darcy flow. 

Darcy flow is a simple model of steady state flow of an incompressible fluid in a porous medium. It 
models the statement that flow is from high to low pressure. For a fixed pressure gradient, the velocity 
is proportional to the permeability k of the medium and inversely proportional to the viscosity jj. of the 
fluid. Let the domain be M, a polygonal planar domain. Assuming that there are no sources of fluid in 
M and there is no other force acting on the fluid, the equations of Darcy flow are 

f + -Vp = and divi^ = inM with v-h = y/ ondM. (11.1) 
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Figure 10: Darcy flow using discrete exterior calculus. The boundary condition is that fluid is coming in from left 
and leaving from right with velocity 1. The velocity inside should be constant and pressure should be linear. The 
flux and pressure are computed in a mixed formulation. The flux is taken to be a primal 1 -cochain associated vnth 
primal edges, and the pressure is a dual 0-cochain on dual vertices, which are circumcenters of the triangles. The 
velocity is obtained by Whitney interpolation of the flux, which is sampled at the barycenters. See Section 11.2 for 
more details. 



where k > is the coefficient of permeability of the medium ;u > is the coefficient of (dynamic) viscosity 
of the fluid, y/ : dM — ► K is the prescribed normal component of the velocity across the boundary, and h 
is the unit outward normal vector to dM. For consistency fgj^y/dT = 0, where dT is the measure on dM. 
Since divograd = A, the simplified Darcy flow equations above are equivalent to Laplace's equation. 

Let ^ be a simplicial complex that triangulates M. Instead of velocity and pressure, we will use flux 
and pressure as the primary unknowns. The flux through the edges is f = * and thus it will be a primal 
1- cochain. Although PyDEC does not implement a flat operator, this is not an issue here because we 
never solve for v, and make / itself one of the unknowns. This implies that the pressure p will be a 
dual 0-cochain since *dp has to be of the same type as /. The choice to put flux on primal edges and 
pressures on circumcenters can be reversed, as shown in a dual formulation in [26] . In exterior calculus 
notation, the PDE in (11.1) is -{jilk) * f + dp = and df = 0, which, when discretized, translates to the 
matrix equation 
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In PyDEC, the construction of this matrix is straightforward. Once a simplicial complex sc has been 
constructed, the following 3 lines construct the matrix in the system above : 

dl = sc[l].d; starl = sc [1] . star 
A = bmat ( [ [ ( -mu/k) * sc [1] . star , sc[l].d.T], 
[sc[l].d, None]], f ormat = ' csr ' ) 

After computing the boundary condition in terms of flux through the boundary edges, the linear 
system is adjusted for the known values and then solved for the fluxes and pressures. Figure 10 shows 
the solution for the case of constant horizontal velocity and linear pressure gradient. 

1 1 .3 Cohomology basis using Hodge decomposition 

The Hodge Decomposition Theorem [1, page 539] states that for a compact boundaryless smooth man- 
ifold M, for any p-form oj e QP{M}, there exists an a e QP"^(M), f5 e n^'^^iM), and a harmonic form 
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h £ D.f (M) such that o) = da + 8(i + h. Here harmonic means that A/i = 0, where A is the Laplace-deRham 
operator d5 + 5d. Moreover da, 5/3 and h are mutually -orthogonal, which makes them uniquely 
determined. In case of a manifold with boundary, the decomposition is similar, with some additional 
boundary conditions. See [1] for details. 

The Hodge-deRham theorem [1], relates the analytical concept of harmonic forms with the topo- 
logical concept of cohomology. For any topological space, the cohomology groups or vector spaces of 
various dimension capture essential topological information about the space [45] . For the manifold M 
above, the p-dimensional cohomology group with real coefficients, which is a finite-dimensional space, 
is denoted HP{M;W) or just H^iM). For example, for a torus, has dimension 2. For a square with 4 
holes used in this example, which does have boundaries, has dimension 4. The elements of HP{M) 
are equivalence classes of closed forms (those whose d is 0). Two closed forms are equivalent if their 
difference is exact (that is, is d of some form). While the representatives of 1 -homology spaces can be 
visualized as loops around holes, handles, and tunnels, those of 1 -cohomology should be visualized as 
fields. If the space of harmonic forms is denoted ^P{M), then the Hodge-deRham theorem says that it 
is isomorphic, as a vector space, to the p-th cohomology space Hf{M) in the case of a closed manifold. 
See [38] for details. Again, the case of M with boundary requires some adjustments in the definitions, as 
given in [1]. 

For finite dimensional spaces, Hodge decomposition follows from very elementary linear algebra. If 
U, V and W are finite-dimensional inner product vector spaces and A:U and B:V are linear 
maps such that BoA = then middle vector space V splits into 3 orthogonal components, which are im A, 
imB^ , and VeiA^ n kerB. In this example, we find a basis for for a square. This is done by finding 
a basis of harmonic 1-cochains. Thus given a 1-cochain a», its discrete Hodge decomposition exists and 
is w = do a -I- ^2 )6 -I- /i. In this example, the cochains a and are obtained by solving the linear systems 
5i do a = 5i w and di 52 )6 = di (o. The harmonic component can then be computed by subtraction. 

In the example code, the main function is the one that computes the Hodge decomposition of a given 
cochain omega. First empty cochains for alpha and beta are created: 

sc = omega . complex 
p = omega . k 

alpha = sc . get_cochain (p - 1) 
beta = sc . get_cochain (p + 1) 

Now the solution for alpha and beta closely follows the above equations for a and j6: 

A = delta (d ( sc . get_cochain_basis (p - l))).v 

b = delta ( omega) . V 

alpha. V = cg( A, b, tol = le-8 ) [0] 

A = d (delta ( sc . get_cochain_basis (p + l))).v 
b = d (omega) . v 

beta.v = cg( A, b, tol=le-8 ) [0] 

Even though the matrices A above are singular, the solutions exist, and since conjugate gradient is used, 
the presence of the nontrivial kernels does not pose any problems [9] . 

The harmonic 1 -forms shown in Figure 11 are obtained by decomposing random 1 -forms and re- 
taining their harmonic components. Since the initial basis has no particular spatial structure, an ad hoc 
orthogonalization procedure is then applied. For each basis vector, the algorithm identifies the compo- 
nent with the maximum magnitude and applies Householder transforms to force the other vectors to 
zero at that same component. 
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Figure 11: Four harmonic cochains form a basis for the first cohomology space H for a mesh with four holes. Each 
cochain is visualized above as a vector field by interpolating it from the edge values using Whitney interpolation. 
See Section 11.3 for details. 

1 1 .4 Sensor network coverage 

As discussed in Section 5, sensor network coverage gaps can be identified with coordinate-free methods 
based on topological properties of the Rips complex. This is for an idealized abstraction of a sensor 
network. The following example constructs a rips_coinplex object from a set of 300 points randomly 
distributed over the unit square, as illustrated in top left in Figure 12. Recall that the Rips complex is 
constructed by adding an edge between each pair each pair of points within a given radius. Top right of 
Figure 12 illustrates the edges of the Rips complex produced by a cut-off radius of 0.15. The triangles of 
the Rips complex, illustrated in bottom left of Figure 12, represent triplets of vertices that form a clique 
in the edge graph of the Rips complex. The Rips complex is created by the following two lines of code 

pts = read_array ( ' SOOpts . mtx ' ) # 300 random points in 2D 
rc = rips_complex ( pts, 0.15 ) 

The sensor network is tested for coverage holes by inspecting the kernel of the matrix Ai = 5f 3i -i- ^2 5j 
[17]. Specifically, null-vectors of Ai, which are called harmonic 1-cochains (by analogy with the defini- 
tion of harmonic cochains used in the previous subsection), reveal the presence of holes in the sensor 
network. In this example we explore the kernel of Ai by generating a random 1-cochain x and extracting 
its harmonic part using a discrete hodge decomposition as outlined in the previous subsection. If the 
harmonic component of x is (numerically) zero then we may conclude with high confidence that Ai is 
nonsingular and that no holes are present. However, in this case the hodge decomposition of x produces 
a nonzero harmonic component h. Indeed, plotting h on the edges of the Rips complex localizes the 
coverage hole, as the bottom right of Figure 12 demonstrates. 

To set up the linear systems, the boundary matrices are obtained from the Rips complex rc created 
above: 

cmplx = rc . chain_coniplex () # boundary operators [ bO , bl , b2 ] 
bl = cmplx [1] . astype (float ) # edge boundary operator 
b2 = cmplx [2] . astype (float ) # face boundary operator 
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Figure 12: Hodge decomposition for finding coverage holes in an idealized sensor network. Top left shows a sam- 
ple sensor network with 300 randomly distributed points. Pairs of points within a fixed distance of one another are 
connected by an edge and this is shown on top right. Triangles are added to the Rips complex when three points 
form a clique (complete graph). These are shown in bottom left. The existence of a harmonic 1-cochain indicates 
a potential hole in the sensor network coverage. In the bottom right figure, edge thickness reflects the magnitude 
of the harmonic cochain on each edge. See Section 11.4 for more details. 



Then the random cochain is created and the Hodge decomposition computed, to find the harmonic 
cochain which is then normalized: 

X = rand (bl . shape [1] ) # random 1-chain 

# Decompose x using discrete Hodge decomposition 
alpha = cg( bl * bl.T, bl * x, tol=le-8)[0] 
beta = cg( b2.T * b2 , b2 . T * x, tol=le-8)[0] 

h = X - (bl.T * alpha) - (b2 * beta) # harmonic component of x 
h /= abs(h).max() # normalize h 



11.5 Least squares ranking on graphs 

This is a formulation for ranking alternatives based on pairwise data. Given is a collection of alternatives 
or objects that have to be ranked, by computing a ranking score that sorts them. The ranking scores 
are to be computed starting from some pairwise comparisons. Some examples of objects to be ranked 
are basketball teams, movies, candidates for a job. Typically, the given data will not have pairwise com- 
parisons for all the possible pairs. There is no geometry in this application, hence no exterior calculus is 
involved. PyDEC however still proves useful because the full version of this example [33] uses an abstract 
simplicial 2-complex. Thus PyDEC is useful in forming the complex and for determining its boundary 
matrix. In this simplified example, only an abstract simplicial 1 -complex is needed. 
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Form a simple graph G, with the objects to be ranked being the vertices and with an edge between 
any two which have pairwise comparison data given. If there are n objects, possibly only a sparse subset 
out of all possible O(n^) pairs may have comparison data associated with them. Here we'll only require 
that the graph be connected. This condition can be dropped with the consequence that the rankings of 
separate components become independent of each other. The comparison values are real numbers. 

Since G is a simple graph, by orienting the edges arbitrarily it becomes an oriented 1 -dimensional 
abstract simplicial complex. The vector a> of pairwise comparison values is a 1-cochain since if A is 
preferred over B by, say, 4 points, then B is preferred over A by -4 points. 

The ranking scores a which are to be computed on vertices form a 0-cochain. For any edge e = 
{u, v) from vertex u to vertex v, the difference of vertex values a{v) - a{u) should match a){e] as much 
as possible, for example, in a least squares sense. This idea is from [40] who proposed it as a method for 
ranking football teams. By including the 3-cliques as triangles, G becomes a 2-dimensional simplicial 
complex. This was used in [37] to extend this ranking idea. In [37] the computation of the scores a is 
interpreted as one part of the Hodge decomposition of to. See Section 11.3 above for a basic discussion 
of Hodge decomposition where it is used for computing harmonic cochains on a mesh. Here we will just 
compute the ranking score a. This is done by solving the least squares problem dl a^co. 

The graph in this example is used for ranking basketball teams, using real data for a small subset of 
American Men's college basketball games from 2010-2011 season. Each team is a node in the graph and 
has been given a number as a name. An edge between two teams indicates that one of more games have 
been played between them. The score difference from these games becomes the input 1-cochain (o, with 
one value on each edge. If multiple games were played by a pair the score differences were added to 
create this data. The data is stored as a matrix in which the first two columns are the teams and the third 
column is the value of the 1-cochain on that edge. 

Once this data is loaded from file, an abstract simplicial complex is created from the first two columns 
which form the edges of the graph. The loading and complex creation is done by the following few lines 
of code 

data = loadtxt ( ' data . txt ' ) • astype ( int ) 
edges = data [ : , : 2] 

# Create abstract simplicial complex from edges 
asc = abstract_simplicial_complex ( [edges] ) 

In PyDEC, the simplices that are given as input to construct a complex are preserved as is. Lower dimen- 
sional simplices that are derived from them are stored and oriented in sorted order. Thus in the above 
data, the edge between node 8 and 1 will be oriented from 8 to 1. The above example data may mean, for 
example, that team labelled 1 lost to team labelled 8 by 9 points. 

The 1-cochain o) is now extracted from the data array and the boundary matrix needed is obtained 
from the complex and the least squares problem solved. All this is accomplished in the following lines 

omega = data[: ,-1] # pairwise comparisons 

Bl = asc . cliaiii_complex () [1] # boundary matrix 

alpha = lsqr(Bl.T, omega) [0] # solve least squares problem 

The resulting alpha values computed are given below. 

Team number 01234567 8 9 10 
a value 14.5 2.0 0.0 7.2 5.9 9.0 2.3 23.6 11.0 23.8 21.7 

The team with a = is the worst team according to these rankings and the one with the largest a value 
(23.8 here) is the best team. Note that the score difference from the game between 8 and 1 happens to 
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Figure 13: A typical graph of a subset of basketball games. Each node is a team, labelled by a number. An edge 
represents one or more games played by the two teams connected by it. The actual graph used in the example in 
Section 11.5 is a complete graph. 

be exactly the difference in a values between them. This won't always be true. See for example teams 3 
and 9. Such discrepancies come from having a residual in the least squares solution, and that is a direct 
result of the presence of cycles in the graph. In fact it may even happen that team A beats B, which beats 
C, which in turn beats A. Thus no assignment of a values will resolve this inconsistency. This is where 
a second least squares problem and hence a Hodge decomposition plays a role. The second problem is 
not considered in this example. See [33] for details on the second least squares problem and the role of 
Hodge decomposition and harmonic cochains in the ranking context. 

12 Conclusions 

PyDEC is intended to be a tool for solving elliptic PDEs formulated in terms of differential forms and 
for exploring computational topology problems. It has been used for numerical experiments for Darcy 
flow [34], computation of harmonic cochains on two and three dimensional meshes [32], and least 
squares ranking on graphs [33]. It has also proved valuable in computational topology work [20, 22], 
for creating complexes and computing boundary matrices. The design goals for PyDEC have been effi- 
ciency and ability to express mathematical formulations easily. Section 1 1 which described some exam- 
ples should give an idea of how close the PyDEC code is to the mathematical formulation of the problems 
considered. Many packages exist and are being created for numerical PDE solutions using differential 
forms. There are also many excellent computational topology packages. PyDEC can handle a large va- 
riety of complexes and provides implementations of discrete exterior calculus and lowest order finite 
element exterior calculus using Whitney forms. These qualities make it a convenient tool to explore the 
interrelationships between topology, geometry, and numerical PDEs. 
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Appendix 

Proof of Propositions.?: By definition Mp(?,j) = Jj^|<W(c7[')*,W((Tj)*>ju. The integrand is nonzero only 

in St(cr^) n St(a^) since the Whitney form corresponding to a simplex is zero outside the star of that 
simplex [21]. Thus 




By the definition of the Whitney map, Mp{i, j) is 
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which completes the proof. 
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